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Abstract 

For axisymmetric models for coronal loops the relationship between the bifur- 
cation points of magnctohydrodynamic (MHD) equilibrium sequences and the 
points of linear ideal MHD instability is investigated imposing line-tied boundary 
conditions. Using a well-studied example based on the Gold-Hoyle equilibrium, 
it is demonstrated that if the equilibrium sequence is calculated using the Grad- 
Shafranov equation, the instability corresponds to the second bifurcation point 
and not the first bifurcation point because the equilibrium boundary conditions 
allow for modes which are excluded from the linear ideal stability analysis. This 
is shown by calculating the bifurcating equilibrium branches and comparing the 
spatial structure of the solutions close to the bifurcation point with the spatial 
structure of the unstable mode. If the equilibrium sequence is calculated using 
Eulcr potentials the first bifurcation point of the Grad-Shafranov case is not 
found, and the first bifurcation point of the Eulcr potential description coincides 
with the ideal instability threshold. An explanation of this results in terms 
of linear bifurcation theory is given and the implications for the use of MHD 
equilibrium bifurcations to explain eruptive phenomena is briefly discussed. 

Keywords: Magnetohydrodynamics; Instabilities; Corona, Structures; Flares, 
Relation to Magnetic Field 

1. Introduction 

Magnctohydrodynamic instabilities of coronal loops are since a long time dis- 
cussed as one of the main theoretical explanations for solar flares, in particular 
compact loop flares (e.g. IPriest[ [1982). Traditionally, investigations of MHD 
instabilities of coronal loops model these loops as straight cylindrical flux tubes of 
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finite length with line tied boundary conditions at the 'photospheric' ends of the 
flux tubes (e.Q. IRaadul[W21IHood and Priest[[W9l[TMTirEinaudi and Van Hovenl | 
19831 IVelli et al\ I1990|) . Such a set-up allows for a wide variety of relatively 
simple equilibrium configurations, hence explaining its popularity. 

The stability of equilibrium configurations of the above mentioned type has 
been studied for several decades using the methods of linear MHD stability anal- 
ysis (e.g. IRaadu| fT972| IHood and Priest] IT9791 fT98T| lEinaudi and Van Hoven] 
H9831 IVelli et al.i [19901 |De Bruyne and Hood[ [T9891 [T9921 IMikic et al.i IT9901 
IHood et aOll994llva"n der Linden and Hoodlll998[[T999|) . In recent years the in- 
vestigations have been extended into the nonlinear regime using large-scale MHD 
simulations (e.ff. |Longbottom et aZ.|H996(|Baty and Heyvaerts[ [l996; Baty[ll99Ta] | 
I1997bl I2000al I2000bl ILionello et al\ ITMHl lArber et al.i \WMi IGerrard et ali 
I200T1 |Browning and Van der Linden[ I2H051 |Browning eta!] 120051 IHood et al.i 
120091) . 

In the present contribution we want to investigate the stability of line-tied 
coronal loop models from a different point of view. The flux tube equilibria 
used to model coronal loops all depend on one or more parameters representing 
quantities like the magnetic twist or the plasma beta. Many investigations study 
how the linear stability of the loops changes as one (or more) of these equilibrium 
parameters vary. 

The systematic variation of one or several parameters of an equilibrium defines 
an equilibrium sequence, and a point of linear instability should correspond to a 
bifurcation point of the equilibrium sequence and vice versa. It has to be kept in 
mind, however, that magnetostatic equilibria are usually calculated by solving 
a mathematically reduced set of equations. It is not at all clear whether there 
is really a one-to-one correspondence between points of linear instability and 
bifurcation points, in particular if line-tied boundary conditions are imposed as 
in models of coronal loops. 

In the present paper we shall investigate the question whether the points of 
linear instability of rotationally symmetric straight line-tied flux tubes have a 
one-to-one correspondence with the bifurcation points of equilibrium sequences. 

We shall use two different ways of calculating the equilibrium sequences, 
namely Grad-Shafranov theory and Eulcr potentials, and wc shall, for simplicity, 
investigate only axisymmetric instabilities and bifurcations. A particularly well- 
studied equilibrium class ( Gold and Hoyle, 1 960) will be used to carry out this 



investigation, mainly because results of linear stability investigations for this 
equilibrium class are readiliy available in the literature (e.g. IHood and Priest! 
IW^fTMTllMikic et aU [TW01|De Bruyne and Hood[ITM2l) . 

In Section [2] the basic equilibrium theory and those parts of the theory of 
linear MHD stability needed in this paper are discussed. The following Section [3] 
presents a brief outline of the numerical method used to calculate the equilibrium 
sequences and to determine their bifurcation points and bifurcating branches. 
The results of these calculations are given in Section 5] and discussed in Section 
[5] The paper closes with a summary in Section [6] 
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2. Basic Theory 

2.1. The Gold-Hoyle Equilibrium 

Wc start our investigation from static equilibrium solutions of the MHD equa- 
tions, i.e. solutions of 

j x B - Vp = (1) 
V x B = Mo j (2) 
VB = 0. (3) 

We are looking for solutions in cylindrical coordinates r, </>, z, and restrict the 
spatial domain to < z < L. The solutions will be considered as straight 
flux tube approximations of coronal loops, in the sense of a large aspect ratio 
expansion. In this case L is the loop length and the boundaries z = and z = L 
have to be identified with the photospheric end points of the loop. The centre of 
the loop is given by r = 0. In the present paper we will only consider solutions 
which do not depend on (j), i.e. axisymmctric solutions. 

We normalise the magnetic field to the value of B z in the centre of the loop 
(r = 0), Bo, the coordinates and the loop length by a typical radial length scale, 
b, and the pressure by Bq/^q. In this normalisation, the Gold-Hoyle equilib- 
rium (Gold and Hoyle, 1960[) is given by the magnetic field components (see e.g. 



|Longbottom et at. I1996|) 



B r = 0, (4) 

B * = 7^2- ( 5 ) 

B, = ~r~ 2 ' (6) 
1 + H 



and the plasma pressure 



1 1-A 2 

The parameter A controls both the field line twist $ between z = and z = L, 

• (8) 

and the plasma beta. For A = 1, the equilibrium is force- free, i.e. the current 
density is parallel to the magnetic field lines, whereas for A = the current 
density is everywhere perpendicular to the magnetic field lines. For values of 
A between and 1 we have a combination of field-aligned and perpendicular 
current density. The equilibrium class is not defined for A > 1 because the 
pressure would become negative in this case. 

The Gold-Hoyle equilibrium class depends only on the variable r and is there- 
fore a one-dimensional MHD equilibrium. One-dimensional equilibria of this type 
can be easily calculated (see e.g. IPriestl 119821 chapter 3.3) We use it here as a 
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kind of prototype flux tube equilibrium, because linear stability results are read- 
ily available for the Gold-Hoy le equilibrium class ( |De Bruyne and Hood, 1992[ ). 
To obtain genuinely two-dimensional equilibria depending on r and z we have 
to resort to one of the more general theories described in the next sections. 

2.2. Grad-Shafranov Theory 

To satisfy the solenoidal condition we write the magnetic field in the form 

B = x e + B^. (9) 

Here the flux function A and the 0-component of the magnetic field depend only 
on r and z. Taking the scalar product of Equation ^ with B and using the fact 
that under the condition of axisymmctry the pressure also depends only on r 
and z, we find that the pressure is a function of A : 

p{r,z)=F(A(r,z)). (10) 

An investigation of the (^-component of |T]) shows that 

bt(r,z)=rB (f ,(r,z,) = G(A(r,z)) (11) 

is also a function of A only. 

The force balance equation can then be reduced to a single partial differential 
equation for A (e.g. IBateman| H978|) : 

^■(\vA)=r^ + hA. (12) 
*■ ' dA r v dA 



r 



The dependence of the pressure and the 0-component of the magnetic field 
on the flux function A have to be specified for a solution of this equation. For 
the Gold-Hoyle solution the flux function A is given by 

A Gii (r,z) = ^\n(l+r 2 ). (13) 

Using equation (|13[) to determine r as a function of A, and substituting this 
expression into Equations ([5]) and ((6]) we find that 

P = i(l-A 2 )cxp(-lA), (14) 

& = 1 -cxp y-^A ■ (15) 

For each value of A, the Gold-Hoyle equilibrium is therefore a solution of the 
Grad-Shafranov equation 

d fldA\ Id 2 A 1 - A 2 / 4 , 

2r — - — exp — —A 



dr \ r dr ) r dz 2 A V A 

2 



1 - exp ( -\a 



A 



-\a\ 
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if the boundary conditions at z = and z = L are given by A = A/2 ln(l + r 2 ) 
as well. As the Gold-Hoylc solutions depend on the parameter A they define 
a solution branch of the Grad-Shafranov Equation (|16[) . As this equation is 
nonlinear it is to be expected that it can also have other solution branches for 
the same boundary conditions. Points where the Gold-Hoyle solution branch 
and any other solution branches meet are called bifurcation points. Standard 



bifurcation theory (e.g. Iooss and Joseph 1980]) tells us that at bifurcation points 



the stability of the solution branches can change. We will discuss this possibility 
in more detail in Section! 

2.3. Euler Potentials 

The Grad-Shafranov theory is very useful for symmetric plasma systems, but a 
Grad-Shafranov type equation can only be derived for translational, rotational 
and helical symmetry (jSolovevl H5571 lEdenstrasserl I1980a"l I1980b|) . Without such 
a symmetry, we have to use a different way to calculate MHD equilibria. The 
approach coming closest to the use of a flux function for symmetric systems is to 
use Euler potentials (sometimes also called Clebsch coordinates) to describe the 
magnetic field. The Euler potential approach has the advantage that it can be 
used to describe three dimensional magnetic fields without symmetry, although 
there are some restrictions concerning the existence of Euler potentials for given 
magnetic fields (e.g. IHessel 119881 IRosner et all I1989|) . For the flux tube like 
equilibria considered in the present paper these constraints do not apply. 

Another reason for using Euler potentials even for symmetric cases is that 
certain types of constraints are a lot easier to impose with Euler potentials 
then with a Grad-Shafranov description. A typical example from solar physics 
is the quasi-static shearing of magnetic arcades. In this case the footpoint dis- 
placement of the fieldlines is the physical parameter which is determined by 
the boundary conditions. In this case the use of Euler potentials is very useful 
(e.g. IBarnes and Sturrockl [W2l |Zwingmann] [15571 IPlatt and Neukirchi IT5M1 
lAntiochos etaL[ HM^ . 

With the Euler potentials a and (3 a general magnetic field can be written as 

B = Vax V/3. (17) 

Any vector field of this form is automatically solcnoidal. In the present paper we 
restrict our analysis to axisymmetric fields. In this case the Euler potential a is a 
function of r and z only, the Euler potential f3 is chosen as /3(r, (f>, z) = (3(r, z) + (f> 
and Equation ([TT]) reduces to 

B = -Va x e<A + Va x V/3. (18) 
r 

By comparison with Equation © , we see that now the Euler potential a corre- 
sponds to the flux function A, whereas has been replaced by the Va x V/3. 
Substitution of Equation ([T5]) into Equation ([1]) and @ gives the two equations 
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for a and (3 {e.g. IZwingmann] [T9871 IPlatt and Neukirchl [T994)) : 



V/3 • V x (Vq x V/3) - V • ( ^Va ) = (19) 

r z I da 



Va ■ V x (V/3 x Va) = 0. (20) 

For the Gold-Hoylc solution «gh is identical with the flux function Aqh given 
in Equation (|13|) . The pressure function p{a) has the same form as p(A) in 
Equation (fTl| . only with a replacing A. We can use to work out that for the 
Gold-Hoyle solution 

/3gh = -\z. (21) 
The function (3 represents the fieldline twist for the Gold-Hoyle solution since 

G n(r,O)-P GIl (r,L) = j = $ (22) 

We impose boundary conditions for both Euler potentials. The boundary 
conditions for a are the same as for A in the Grad-Shafranov case. For (3 we 
use Equation ([2"T]) on the boundaries. This fixes the footpoint displacement of 
ficldlincs crossing the boundaries z — and z = L. In the same way as in 
the Grad-Shafranov case the Gold-Hoyle solutions arc a solution branch for the 
Equations (HHJ) and (|2T)|) . The same statements about bifurcations and stability 
apply as in the Grad-Shafranov case. 

2.4. Linear Stability 

The theory of linear MHD stability is a vast area and we only summarise some 
results which are important for the following discussion. Defining the Lagrangian 
displacement £ the linearized ideal MHD equations can be written as 

Po^| =F«) (23) 

where 

F(€) = — [(V x Br) x B ] + — [(V x B ) x B x ] + V(£ • Vp +7PoV • £)• (24) 
Mo Mo 

The components of the magnetic field perturbation Bi are given by 

B 1= Vx(£xB ). (25) 

For the problem we are discussing in the present paper Equation (f2"3"]) has to 
be solved on a tube-like domain with line-tying boundary conditions at z = 
and z = L, with Bo and po given by the Gold-Hoyle solution. The line-tying 
condition corresponds to 

£ = (26) 
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on the boundaries. Assuming an exponential time-dependence for the pertur- 
bation one obtains a self-adjoined eigenvalue problem. Instabilities occur 
when one of the eigenvalues of the equation changes sign. The corresponding 
perturbations £ can be classified according to their different spatial structure. 
In general we speak of different modes when refering to the spatial structure of 
the instabilities. For each mode it is usually sufficient to investigate the largest 
eigenvalue corresponding to this mode, as it is this eigenvalue which determines 
whether a mode is stable or unstable. 

It can be shown that the Gold-Hoyle solution is always stable for A = 1 , i.e. in 
the force-free case. When decreasing A the value of A where the Gold-Hoyle solu- 
tions become unstable to the different possible modes under line-tying boundary 
conditions depends on the length of the flux tube L. A thorough investigation of 
this problem has been carried out by |De Bruyne and Hood| (|1992p and we will 
make use of their results in the later parts of this paper. Since we investigate 
only axisymmctric equilibria and bifurcations we will also restrict our attention 
to the axisymmetric modes (sometimes called "sausage modes" ) . 



3. Numerical Method 

The numerical calculations have been carried out with a code based on a continu- 
ation method (e.g. |Allgower and Gc org, 1990J). The code used here is based on a 
method proposed by(Keller ( 197TJ) and has been successfully applied to a variety 
of problems in plasma physics, solar physics, magnetospheric physics and astro- 

] [Tg53inWllNeukirchlll993aU1993birNeukirch and Hessel | 



physics (e.g. Zwingmann 



iPlat t and Neukirc 



El [1194; Schr oer et al\ IT9M1 IBecker et al\ [TMcIl IMJTl 



IRomeou and Neukirch[[T9M[M?Tl[M^all2002bl|Kiessling and Neukirch[[M?3l) . 

The method has the advantage that it can calculate sequences of equilibria 
depending on an external parameter (like A for the Gold-Hoyle solutions), and 
detect bifurcation points. It is also possible to calculate bifurcating equilibrium 
sequences. The code uses a finite element discretization allowing for a flexible 
grid structure. Further details can be found in INeukirchl (|1993a[) and INeukirchl 
(|1993bl) . 

We have solved both the Grad-Shafranov Equation (16]) and the Euler po- 
tential Equations (|19p and (|20[) on a numerical domain extending from r = to 
r = 8 and from z = to z = L, where L is varied between 3 and 8. The radial 
extent of the domain is chosen along the same lines as done by |Longbottom et al. 
(|I996p in their MHD simulations of the sausage instability. 

For the Grad-Shafranov equilibrium sequences we have used 

A b = A GK (27) 

as boundary condition on all boundaries. In the Euler potential case the bound- 
ary conditions are given by 

a b = A GU , (28) 
fa = /3gh- (29) 
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Note that in both cases both the differential equation and the boundary condi- 
tions depend on the parameter A. Due to the boundary conditions the Gold-Hoyle 
solutions are one solution branch of the equations. This can be used to check the 
accuracy of the numerical code and to adjust the resolution. In all runs presented 
in Section |4] we have used a numerical grid with 1800 triangular finite elements 
corresponding to a resolution of 61 by 61 grid points in each spatial direction. 
The grid is equidistant in the z-direction but non-equidistant in the r-direction 
with a higher resolution towards the axis of the tube. 



4. Results 

For both the Grad-Shafranov and the Euler potential case we have carried 
out a numerical investigation of the bifurcation properties of the Gold-Hoyle 
solution branch using the numerical method described in Section [3l For a series 
of values of the loop length L, we have first calculated the Gold-Hoyle branch 
with our code, starting with the force- free solution (A = 1) and then following 
the branch for decreasing A into the non-force- free regime. Although we know the 
Gold-Hoyle branch analytically this procedure allows us to check the accuracy 
of our numerical calculations and to use the capability of the code to detect 
bifurcation points. At such points other solution branches cross the Gold-Hoyle 
branch, and we expect that those points correspond to the instability threshold 
of the m = 0-instability under line-tying conditions as, for example, calculated 
by |De Bruyne and Hood (|1992[) . For the detected bifurcation points we have 



then also calculated the bifurcating branches for a range of A values. This is 
important to check whether the spatial structure of the bifurcating solution 
branch coincides with the predictions made by linear stability theory on the 
basis of the structure of the unstable mode. 

An important point to emphasize here is that because we calculate the Gold- 
Hoyle branch in the direction of decreasing A, we will also number the bifurcation 
points in this direction, i.e. when we speak of first and second bifurcation the 
A value of the first bifurcation will be bigger than the A value of the sec- 
ond bifurcation. Although this is opposite to the terminology normally used 
in bifurcation theory, we have decided to keep the parametrization used by 
De Bruyne and Hood| (|1992| to make a comparison with their results easier. 



4.1. The Grad-Shafranov Case 

We have carried out calculations of the Gold-Hoyle branch for loop lenghts L = 
3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 and 10.0. In all cases, we have calculated the 
Gold-Hoyle branch until we had found at least two bifurcation points. The A 
values of the two first bifurcation points found by code for the different L values 
are listed in Table [TJ A graphical representation of these values is shown in 
Figure [1] In this figure we plot the values for the first and second bifurcation 
points in the L-A-plane. Also shown in Figure Q] is the stability threshold for the 
0-instability (dashed line) derived by |De Bruyne and Hood (|1992p using 



linear MHD theory. The figure clearly shows that the linear stability threshold 
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Table 1. The first and second bifurcation points for the Grad-Shafranov 



case. 


L 


3.0 


4.0 


5.0 


6.0 


7.0 


8.0 


9.0 


10.0 


Ai 
A 2 


0.598 
0.472 


0.629 
0.528 


0.652 
0.566 


0.663 
0.593 


0.675 
0.612 


0.684 
0.627 


0.684 
0.639 


0.690 
0.648 



1 .0 
0.8 

0.6 

0.4 

0.2 
0.0 

2 4 6 8 10 

L 



Figure 1. The dependence of the A values of the first (0) and the second (x) bifurcation point 
on the loop length L for the Grad-Shafranov case. The dashed line is the instability threshold 
of the m = mode derived by [De Bruync and Hood (1992 ) using linear MHD stability theory 
under line-tying conditions. Solutions with A values below the dashed line are unstable with 
respect to the sausage mode. It is obvious that the second and not the first bifurcation point 
for a given loop length corresponds to the m = instability. 

corresponds to the second bifurcation along the Gold-Hoylc branch. This raises 
the question what the first bifurcation point corresponds to. 

To answer this question, we have calculated the bifurcating branches for the 
first and second bifurcation points for loop lengths of L = 3.0, 5.0 and 7.0. The 
structure of the bifurcation diagrams is very similar for all three cases and we 
therefore only show the case L = 7.0 (Figure^. The four quantities shown in 
Figure [Hare the polodial magnetic energy 

Wp = / \ (} VA ) dV > (30) 



SOLA: tnromeou_revised.tex; 30 October 2009; 14:56; p. 9 



© o 

w — — X 



© 
x " 



© 
x - 



© 
•x 



-X- 



Neukirch and Romeou 



2.0 




Figure 2. Bifurcation diagrams of the Grad-Shafranov case for L = 7.0. Shown are the 
poloidal magnetic energy (upper left), the toroidal magnetic energy (upper right), the thermal 
energy (lower left), and the free energy as defined by Grad ( 1964 ) (lower right). For definitions 
of these quantities see the main text. The free energy is only shown for values of A close to the 
second bifurcation point to make the difference between the branches more obvious. 

the toroidal magnetic energy 

Wt = Jl (h<\ dV, (31) 

the thermal energy 

Wt h = J P dV, (32) 
and the free energy defined by IGradl (|1964p 

W f = W p - (W t +W th ). (33) 

At the first bifurcation point another solution branch crosses the Gold-Hoyle 
branch. The bifurcating branch exists for values of A both smaller and larger 
than the bifurcation point A. At the bifurcation point the poloidal and toroidal 
magnetic energies of the bifurcating branch go from values smaller than the 
energies of the Gold-Hoyle branch to values larger than the Gold-Hoyle branch 
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A = 0.646457 A = 0.600544 



2 4 6 8 




Figure 3. Left: Solution for A = 0.646457 on first bifurcating branch. Right: Solution for 
A = 0.600544 on second bifurcating branch. Whereas the solutions on the first bifurcating 
branch show a sin(7r,z/L)-dependence superimposed on the Gold-Hoyle solution, the solutions 
on the second bifurcating branch have a sin(27r,z/L)-dependence. 



in the direction of decreasing A. The thermal energy of the bifurcating branch 
is higher than that of the Gold-Hoyle branch for A larger than the bifurcation A 
and smaller than the Gold-Hoyle branch beyond the bifurcation point. 

The second bifurcating solution branch only exists for values of A which are 
smaller than the bifurcation A. This is to be expected on the basis of standard 
bifurcation theory (e.g. |Iooss and Joseph] [TOSOp . taking the spatial structure of 
the solutions along the branch into account (see below) . The mathematical argu- 
mentation is given in the Appendix. The poloidal and toroidal magnetic energies 
along this branch are larger than those of the Gold-Hoyle branch, whereas the 
thermal energy is smaller than the thermal energy of the Gold-Hoyle branch for 
the same value of A. 

Also shown in Figure is a plot of the free energy for values of A close to 
the second bifurcation point to enhance the difference between the branches. We 
can see that in this range the Gold-Hoyle branch has the biggest free energy. 
The first bifurcating branch has lower free energy than the second branch but 
both branches have a lower free energy than the Gold-Hoyle branch. This shows 
that a transition to the bifurcating branches at fixed A is indeed energetically 
favourable for the system, as the system will always try to settle into a state of 
lower free energy. 

The spatial structure of the solutions on the bifurcating branches is shown in 
Figure [3J The obvious difference between the solutions on the two branches is 
their dependence on z. Whereas the solutions on the first bifurcating branch show 
a sin(7rz/L)-dependence superimposed on the Gold-Hoyle solution, the solutions 
on the second bifurcating branch have a sin(27rz/L)-dependence. Both functions 
are consistent with the boundary condition A = Aqh at z = and z = L. We 
will discuss the implications of this finding in the light of linear stability theory 
in Section [5] 
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Table 2. The first bifurcation point for the Euler potential case. 



L 


3.0 


4.0 


5.0 


6.0 


7.0 


8.0 


9.0 


10.0 


Ai 


0.472 


0.528 


0.566 


0.593 


0.613 


0.628 


0.639 


0.648 



1 .0 
0.8 

0.6 

0.4 

0.2 
0.0 



1 



Figure 4. The dependence of the A values on the loop length L for the Euler potential case. 
The Euler potential case is different from the Grad-Shafranov case because in this case the 
first (0) bifurcation point corresponds to the point where the m = 0-mode becomes unstable. 



4.2. The Euler Potential Case 

In the Euler potential case we have carried calculations of the Gold-Hoylc branch 
for the same values L as in the Grad-Shafranov case. The calculations were run 
for about the same A range as for the Grad-Shafranov equation, but only one 
bifurcation point was detected in this range. The A values of the bifurcation 
point for all loop lengths is given in Table O A comparison with Table [T] shows 
that the first bifurcation point in the Euler potential case corresponds to the 
second bifurcation point of the Grad-Shafranov case. 

This is also obvious if we plot the A values of the bifurcation point for different 
L in the L-A-plane to compare with the results of |De Bruyne and Hood| (j!992p 
(see Figure |4}. It can clearly be seen that in the Euler potential case it is 
obviously the first bifurcation which coincides with the stability threshold of 
linear MHD. This difference between the Grad-Shafranov case and the Euler 
potential case is surprising and we will discuss the reasons for this in Section [5l 
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Figure 5. Bifurcation diagrams of the Eulcr potential case for L = 7.0. Shown are the poloidal 
magnetic energy (upper left), the toroidal magnetic energy (upper right), the thermal energy 
(lower left) and free energy (lower right). The major differences to Figure [2] are that the 
bifurcation point shown in this diagram corresponds to the second bifurcation point of the 
Grad-Shafranov case, and that the bifurcating solution sequence branches off in the direction 
of increasing A. This opposite to the Grad-Shafranov case. 



Another major difference between the Grad-Shafranov case and the Eulcr 
potential case is that in the Euler potential case the new solution sequence 
branches off towards increasing values of A, whereas for the Grad-Shafranov 
case the bifurcating sequence branching off towards decreasing A values. This 
has implication for the stability of the bifurcating branch. 

For the Euler potential case we define the poloidal magnetic energy as 

W p = J\ (} Va ) (34) 
the toroidal magnetic energy as 

W t = J - (Va x V^) dV, (35) 
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Figure 6. Contour plot of the Euler potentials a (left) and (3 (right) on the bifurcating 
branch for A = 0.657. The spatial structure of a is similar to the spatial structure of A in the 
Grad-Shafranov case. 



the thermal energy as 

W th = J pdV, (36) 

and the free energy as 

Wi = W p + W t - W th . (37) 

One should note that the contribution of the toroidal magnetic energy to the 
free energy is positive for the Euler potential case, whereas it is negative in the 
Grad-Shafranov case. The reason for this are the different constraints on the 
system in the two cases (see[Grad ( 1964|). The poloidal magnetic energy of the 
bifurcating branch is slightly lower than that of the Gold-Hoyle solution, whereas 
the toroidal magnetic energy is higher. The thermal energy is also higher than 
the thermal energy of the Gold-Hoyle solution and the free energy is slightly 
larger than that of the Gold-Hoyle solution. 

The spatial structure of a on the bifurcating branch is similar to the spatial 
structure of A in the Grad-Shafranov case (see Figure O obviously there is 
no analogue for (3 in the Grad-Shafranov case). In the linear regime a has a 
sin(27rz/L) z-dependence, whereas (3 has a 1 — cos(27r z/L) structure due to 
Equation ([20]) . 



5. Discussion 



The results presented in Section 3] raise the following questions. 

• Why does the m = 0-instability correspond to the second and not to the 
first bifurcation point in the Grad Shafranove case ? 

• Why is the Euler potential case different from the Grad-Shafranov case ? 
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To answer these questions we first analyse the connection between the line- 
tying condition in linear ideal MHD stability and the boundary conditions in the 
Grad-Shafranov and Eulcr potential cases. Close to the bifurcation points, wc 
can represent the solutions on the bifurcating branches for the Grad-Shafranov 
case by 

Aba (r, z, A) = A GH (r, A) + eAi(r, z) + . . . , (38) 

where e <C 1. Note that this expansion differs slightly from the expansion used 
in the Appendix. To first order in e the function A\ has to satisfy the equation 




+ 



i d*& 



2 h 2 



a gh 2r dA 2 



Age 



Ai, (39) 



where p and are given by Equations (| 14[) and (|15[) . Since the boundary con- 
ditions for Abif are already satisfied by Aqu, the function Ai must vanish on all 
boundaries. Since all coefficients of Equation (|39[) depend only on r, it is easy 
to see that the solutions of Equation (f3"9"|) must be have the form 

Ai(r,z) = F n (r) sin (nnz/L), n=l,2,3,..., (40) 

with F n (r) a radial function. 

Close to the bifurcation point, a linear stability analysis of the fundamental 
Gold-Hoyle branch would give Lagrangian perturbations £ which are related to 
Ax through Equation Since the deviation of the poloidal field from the 

Gold-Hoyle branch along the bifurcating branch is given by V x (AiV0), one 
can easily see that 

AiV0 = -(£ • VA GH )W> (41) 

so that 

M = (42) 



since Agh depends only on r. Equation (|42[) shows that the boundary condition 
Ai = only implies £ r = 0, but not necessarily £ = 0. We therefore surmise 
that the first bifurcation in the Grad-Shafranov case corresponds to a linear 
Lagrangian displacement with non- vanishing £^ and/or £ 2 . The solutions on the 
first branch would satisfy the boundary condition A\ = 0, but not £ = 0. The 
first bifurcating branch therefore corresponds to solutions which do not satisfy 
the line-tying boundary conditions, and, for example, can only be reached from 
the Gold-Hoyle branch if flow through the boundary is allowed (£ z ^ 0). 

The second bifurcating branch satisfies both A\ = and £ = 0. This is 
corroborated by the fact that the z-dependence of the r-component of the 
Lagrangian perturbation for the sausage mode is given by sin(27rz/L) {e.g. 



Longbottom et ai\ I1996[) . This matches exactly the z-depcndcncc of A\ on the 



second bifurcating branch. Therefore the bifurcation points and the linear insta- 
bility thresholds coincide. 

In the Euler potential case, we have boundary conditions for both a and (3, 
thus constraining the system more than in the Grad-Shafranov case. One can 
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derive the connection between the Lagrangian perturbation £ and the linear 
perturbations cti and /3i of the Euler potential from the expression for the linear 
perturbation of the magnetic field, 

Bi = V x (£ x B ) = Vat x V0 + V«i x V/3 + Va x V/?i. (43) 

With a bit of algebra one can show (e.g. [Zwingmann. 1987) that 

ai = -£-Va , A = -I- V^ + A,), (44) 

which for axisymmctric flux tube equilibria like the Gold-Hoylc equilibrium 
discussed in this paper leads to 

As is to be expected the expression connecting a\ and £ r is the same as for A\ 
and £ r in the Grad-Shafranov case, and a\ = on the boundaries ensures that 
£ r = on the boundaries. The boundary condition (3\ = imposes an additional 
constraint, which links £^ and £ z on the boundaries, ensuring that £j_ vanishes 
on the boundaries. This is consistent with the line-tying boundary conditions im- 
posed by |De Bruyne and Hood| (|1992|) and explains why the bifurcation points 
coincide with the linear stability threshold in the Euler potential case. 

We suspect, but cannot prove, that the different structure of the bifurcation 
diagrams present in Figures [2] and [5] is also due to the different constraints 
imposed upon the system by using different descriptions for the magnetic field. 
The different structure of the bifurcation diagrams may have implication for the 
stability of the bifurcating equilibrium branch. Usually, when moving along a 
stable equilibrium sequence and crossing a bifurcation point so that the equilib- 
rium sequence is unstable beyond the bifurcation point, the bifurcating branch 
is linearly stable close to the bifurcation point if it bifurcates in the forward 
direction and linearly unstable if it bifurcates in the backward direction (see 
e.g. Iooss and Joseph I1980|) . In the present case this would imply that in the 
Grad-Shafranov case the second bifurcating branch is linearly stable, whereas 
this branch is unstable in the Euler potential case. This is also supported by the 
fact that the second bifurcation branch has a lower free energy than the Gold- 
Hoylc branch for the Grad-Shafranov case, whereas it has a higher free-energy in 
the Euler potential case. It has to be remarked, however, that this is a conjecture 
as we have no rigorous proof. 



6. Summary and Conclusion 

We have investigated the relationship between MHD bifurcation and linear sta- 
bility for a class of axisymmctric straight flux tubes under line-tying boundary 
conditions. For simplicity we only considered rotationally symmetric pertur- 
bations, allowing only for sausage modes. We have used two different ways 
of calculating the equilibrium sequences including bifurcating branches - one 
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approach uses the Grad-Shafranov equation, the other approach uses Euler 
potentials. It turns out that only the Euler potential case shows a one-to-one 
correspondence between the first bifurcation point and the linear instability 
threshold for the sausage mode. The Grad-Shafranov case shows an additional bi- 
furcation which does not correspond to the instability threshold under line-tying 
boundary conditions. This difference can be explained by the different con- 
straints imposed on the bifurcating equilibrium branches in the Grad-Shafranov 
and the Euler potential cases. 

Furthermore, even though the second bifurcation point of the Grad-Shafranov 
case coincides with the first bifurcation point of the Euler potential case and 
the linear instability threshold, the structure of the bifurcation diagrams differ 
considerably between the Grad-Shafranov and the Euler potential case. The 
reason for this is not yet clear, but is probably also due to the difference in 
boundary conditions. In any case this difference has implications for the stability 
of the bifurcating equilibrium branches (see e.g. |Iooss and Jose ph, 19cTD| and is 
therefore important to decide whether the system is able to find a new equilib- 
rium (in the present case a new axisymmctric equilibrium) if one would consider 
an imaginary process driving the flux tube across the instability threshold. 

The present investigation is a preparation for studying equilibrium sequences 
of magnetic flux tubes and other solar magnetic structures together with their 
bifurcations in three dimensions. Preliminary steps have already been made (see 
e.q. lRomeou and Neukirchl I2002a|) and more detailed investigations are planned 
for the future. 
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Appendix 

Whereas the first bifurcating branch in the Grad-Shafranov case exists for values 
of A which are both bigger and smaller than the A at the bifurcation point, 
the second bifurcating branch exists only for A smaller than the bifurcation A. 
This fact can be explained by using standard bifurcation theory to calculate 
the structure of the bifurcating branches close to the bifurcation points. The 
argument is actually independent of the form of the functions p(A) and b^(A). 
The qualitative structure of the bifurcation diagram will thus be the same even 
if p(A) and b^(A) are changed as long as the fundamental branch consists of 
solutions which depend only on the radial coordinate r. 

We start by writing the Grad-Shafranov equation in the form 




(46) 
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where the function N(A, r, A) summarizes the nonlinear part of the Grad-Shafra- 
nov equation given by p(A,X) and b$(A,X). For the present paper p and 
are given by Equations (|14[) and (|15|) . For the following argument, however, 
the exact form of N(A, r, A) is irrelevant, as long as it is analytic in A and 
A at the bifurcation points we want to investigate. We will not give here any 
details of the mathematical background which can be found for example in 
IHesse and Schindl er (1986} and |Hesse and K iessling (1987) . These papers treat 
slighly different bifurcation problems, but we will be using the same technique. 

Let A* be the value of A at either of the bifurcation points and let Aq = 
Aq(X*) be the solution of Equation (|46|) at the bifurcation point. To calculate 
the bifurcating branch we expand A and A as 



A = J2e k X k , (47) 

fc=0 
oo 

A = 5> fc A fc , (48) 



fe=0 



where Ao = A* and An as above. Since G is analytic in both A and A for A > 
we can expand Equation (|46|) in a power series in e: 

° = E^W)' A (^ r ) e= /- ( 49 ) 



fc! de k 

k=0 



As each power of e must satisfy this equation independently we obtain 

i-G(A(e),A(e),r) =0, k = 0, 1, 2, . . . . (50) 



de 



£=0 



Obviously, the lowest order equation 

G(A (\*),\*,r)=0 (51) 

is just the Grad-Shafranov equation at the bifurcation point and therefore triv- 
ially satisfied. 

For the discussion of the higher order equations we first have to look at 
the boundary conditions the A/~ have to satisfy. The boundary condtion Ab for 
A(r 1 z,X) is given by the fundamental branch solution Ab(r, z, A) = Ao(r,z,X) 
(the Gold- Hoy le solution in the present paper). Therefore we can extend Ab 
into the domain. Using the expansion (|4T|) in Ab(r, z, A(e)) we can see that the 
boundary condition each of the A k in Equation l]48p has to satisfy is given by 

A{ » k) = h^ Ao{r > z > x{£)) (52) 



£ = 



Note that A b satifies the same Equation (|50|) as A k - 
For 0(e) we get from Equation (|50| 



G yl (A (A*),A*,r)^ 1 +G A (^ (A*),A*,r)A 1 =0, (53) 
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with 

/ 1 \ dN 
G A A l = -rV ■ (-VA,) - —(A ,X*,r)A 1 , (54) 

dN 

Gx = —(A ,\*,r). (55) 

As mentioned above A]/ satisfies the same equation as A\ and therefore the 
funct 



ion 



A[ =A X - Al 11 (56) 



satisfies GaA[ = or explicitely 



r 



with = on the boundaries. Since all coefficients of Equation ([57)) depend 
only on r its solution can be obtained by separation of variables with the general 
form of A[ being 

A[(r,z) = F n (r)sm(mrz/L), n = 1, 2, 3, . . . . (58) 

The exact form of F n (r) is of no importance for the following argument. 

If we want to calculate Ai, we have to go to the next order (0(e 2 )) of the 
expansion, giving 

„ ( 1 „ , \ dN . ld 2 N „, d 2 N . x 
~ rV ■ J ~ ~dA 2 _ 2^ - 9l5A^ lAl 

l<9 2 iV , dN 

-2lu^-a^ = < 59 > 

where all derivatives of N(A, A, r) are evaluated at the bifurcation point (e = 0). 
Similarly to A^' at O(e), A?^ satisfies the same equation as A 2 . We define 

A' 2 = A 2 - A^ (60) 

which obeys the equation 

GaA< 2 = JSk 2 + 2A^A[) + X^A[. (61) 



By Fredholm's alternative the right hand side of Equation ((5T|) has to be orthog- 
onal to A[ , i. e. 

jf I"" GS'- 4 " + 2A » 1>a; > + A 'lsH • 4 '"* fa " °- (62) 

To proceed we assume in agreement with the Gold-Hoyle solution that the 
function A^ has the form 

4 X) = ^h{r) (63) 
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where fb(r) is left unspecified here. Equation (f62|) can then be used to calculate 
Ai in the form 



Xl l Jo {dA^ fb{r) + dAdx) A[2rdrdz = 

r L r rma "" ld 2 N q 

The double integral on the right hand side of Equation (|64|) can be split into 
two separate integrations over r and z, since the integrand depends on z only 
through A\. As A[ has the form ([58]) . the integral over z is given by 

l L L 
sin 3 (nnz/L)dz = (cosn7r — 1) + - — (cos 3 n7r — 1), (65) 

717T 3?17T 





which vanishes for all even n. As the integral on the left hand side of Equation 
(|64[) is nonzero, this implies that for even n (and in particular for n = 2) Ai 
vanishes. The bifurcation at bifurcation points with modes having even n is 
therefore quadratic. 

This explains the structure of the bifurcation diagram in Figure [21 because 
the first bifurcation obviously corresponds to the A\ for n = 1, whereas the 
second bifurcation corresponds to n = 2. Therefore the structure of the first 
branch close to the bifurcation point is given by 

A = A*+eAi + ..., (66) 
A = A (r,z,X*) + eA 1 (r,z,X*) + .... (67) 

The slope of the bifurcating branch at the bifurcation point is determined by 
Ai =/= in this case and it is obvious that the bifurcating branch exists for both 
A > A* (eAi > 0) and A < A* (eA x < 0). 

Close to the second bifurcation point we have 

A = A* + ie 2 A 2 + ..., (68) 
A = A (r,z,X*)+eA 1 (r,z,X*) + (69) 

because here Ai vanishes. As the correction to A* depends quadratically on e, 
positive and negative e give the same value of A. This implies that the second 
bifurcating branch actually consists of two branches, one for positive and one 
for negative e. Since a change of sign of e in Equation ([69]) corresponds to a 
simple mirroring of the sin(27rz/L) function at the point z = L/2, the two 
branches have exactly the same energies. We remark that since Ai = in this 
case A\ = A\ as the boundary contribution to A[ vanishes. The numerical 
calculations corroborate these results as the same second bifurcation branch is 
found by the code starting both with negative and positive e. The only difference 
between the calculations is the mirroring of the z-dcpendcncc of A along the 
bifurcating branch. 
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